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Gravitational wave sources are a promising cosmological standard candle because their intrinsic 
luminosities are determined by fundamental physics (and are insensitive to dust extinction). They 
are, however, affected by weak lensing magnification due to the gravitational lensing from structures 
along the line of sight. This lensing is a source of uncertainty in the distance determination, even 
in the limit of perfect standard candle measurements. It is commonly believed that the uncertainty 
in the distance to an ensemble of gravitational wave sources is limited by the standard deviation of 
the lensing magnification distribution divided by the square root of the number of sources. Here we 
show that by exploiting the non-Gaussian nature of the lensing magnification distribution, we can 
improve this distance determination, typically by a factor of 2-3; we provide a fitting formula for the 
effective distance accuracy as a function of redshift for sources where the lensing noise dominates. 

PACS numbers: 04.30.Tv,95.36.+x,98.62.Sb,98.80.Es 



I. INTRODUCTION 

In 1998 two groups measuring the luminosity distance- 
redshift relation (the Hubble diagram) from Type la su- 
pernovae (SNe) reported that the expansion of the Uni- 
verse was accelerating [l|, 0] ■ This discovery has stimu- 
lated a range of efforts to measure the cosmic expansion 
history and assess whether it is consistent with a cos- 
mological constant or if alternatives such as quintessence 
are required. The Type la SNe continue to provide one 
of the most valuable constraints Q, due to quality data 
at a range of redshifts, and the lack of cosmic variance 
limitations that plague alternatives such as weak grav- 
itational lensing (WL) and baryon-acoustic oscillations 
(BAO) at low redshift. 

The advent of gravitational wave astronomy has 
prompted interest in gravitational wave (GW) sources as 
a standard candle. Schutz 0] showed that a gravitational 
waveform from merging compact objects can be used to 
measure the distance to the source; a redshift obtained 
of an electromagnetic counterpart or host galaxy would 
then allow one to place the GW source on a luminosity 
distance-redshift relation. The GW source method has 
the key advantage over other standard candles that its 
luminosity can be determined from fundamental physics, 
thus alleviating the common concern with standard can- 
dles that they could evolve with cosmic time. They 
are also insensitive to dust opacity. Finally, many pro- 
posed space-based gravitational wave detectors measure 
test mass separations directly in units of the laser wave- 
length, as opposed to supernovae, which measure rel- 
ative distances and require an independent calibration 
ladder in order to measure absolute distances. Thus GW 
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sources offer the potential for a low-systcmatics probe of 
the expansion history of the Universe. Examples of such 
sources would include mergers of massive black holes, 
observable with the Laser Interferometer Space Antenna 
(LISA [41]), and binary neutron star or stellar mass black 
hole binaries, observable with a more futuristic space- 
based detector in the ~ 1 Hz band such as the Big Bang 
Observer (BBO). The latter, in particular, could poten- 
tially observe hundreds of thousands of neutron star- 
neutron star binary inspirals [5] . 

A final advantage is that the distance determination 
to a GW source is limited by the signal-to-noise of the 
measurement (and by partial degeneracies with binary 
parameters), as opposed to Type la SNe, which have 
a seemingly random ~ 15% scatter in their luminosi- 
ties even after the light curve stretch correction. This 
means that the statistical power of GW sources may 
be interesting even if the number of usable sources is 
far less than the number of usable SNe. In fact, for 
high signal-to-noise detections of GW sources, the dis- 
tance determination should be limited not by the intrin- 
sic dispersion of the source nor by the measurement er- 
ror, but rather by weak lensing magnification: the inter- 
vening matter between us and the source will magnify 
the GW source and affect our measurement of the dis- 
tance. The apparent flux from the source is increased 
by some factor /i, which is usually ~ 1, and the ap- 
parent distance Z? app thus differs from the true distance 
D according to D app — D^ 1 / 2 . This phenomenon of 
course occurs for all standard candles, and has long been 
recognized as an issue for SNe but its importance 

relative to intrinsic dispersion is much greater for gravi- 
tational waves. (Gravitational wave measurements with 
nearby sources or with lower signal-to-noise per source, 
such as the proposed binary progenitors of short gamma 
ray bursts [9|, [l(| , are much less affected since the lensing 
scatter is subdominant.) 
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The usual way of accounting for the magnification ef- 
fect is to suppose that it adds in quadrature with the 
intrinsic-luminosity and apparent flux measurement con- 
tributions to the distance error. That is, 
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where D is the distance, L is the intrinsic luminosity, F 
is the measured flux, and [i is the magnification; with TV 
sources, this uncertainty is reduced by a factor of \jN 
[a El ■ Since the last term dominates for GW sources 
and is significant for high-z SNe, there is great motiva- 
tion to try to reduce it. One possibility would be to try 
to construct an estimated magnification fi from exter- 
nal data sets; the last term should then be replaced by 
a 2 (/i — fi). Unfortunately, WL shear maps have too much 
high spatial frequency noise to be useful as a magnifica- 
tion estimator for point sources [12j , but galaxy maps are 
highly correlated with the mass distribution and may be 
able to reduce the lensing dispersion term by a factor of 
~ \/3 [l3l - [l~5j ]. Maps of flexion (i.e. the gradient of the 
shear measured using the banana- or trefoil-shaped dis- 
tortion of a galaxy 16]) could also be useful if very hig h 
source densities (> lOOarcmin -2 ) can be obtained [17] . 

The purpose of this paper is to explore yet another 
method of reducing the lensing dispersion in Eq. ([1}. 
Because the probability density function (PDF) of fj, is 
highly non-Gaussian (technically non-r, as we show in 
Appendix [Aj , our ability to centroid it using N sources 
can be much better than a^/y/W. This is in fact a famil- 
iar result: as an extreme case, many distributions such 
as the Airy diffraction pattern in optical astronomy can 
be centroided even though their variances are formally 
infinite. In the case of lensing magnification, the vari- 
ance of /x is often dominated by a long tail to positive 
values, corresponding to lines of sight that pass through 
a galaxy or its halo, whereas the information content is 
dominated by relatively empty lines of sight with \i — 1 
sharply peaked around a slightly negative value. In such 
cases, the use of outlier-rejecting statistics not only re- 
moves sources with misidentified redshifts Q but also 
reduces the uncertainty in D(z) for correctly identified 
sources. (A similar point has been made in the recent 
paper by Shang & Haiman fig.) 

This paper is organized as follows. In Section [HI we 
discuss the information-theoretic limits on centroiding a 
distribution. Numerical results and simulations are pre- 
sented in Section Mil Section IIVI gives cosmological con- 
straints obtainable with the reduced centroid errors for 
BBO and for LISA. We conclude and briefly discuss sys- 
tematics in Section IVl 

We focus here on the problem of measuring D(z) from 
GW sources. However, the formalism is applicable to any 
standard candle, and we briefly discuss the implications 
for Type la supernovae. 



II. CENTROIDING A DISTRIBUTION 

In this section, we consider the problem of determin- 
ing the distance D to a population of A" standard sources 
at some redshift z. For simplicity, we consider first the 
case with no intrinsic dispersion in the source luminosity, 
and then generalize to the case with a known additional 
source of scatter (e.g. an intrinsic dispersion or measure- 
ment uncertainty). For large N, the maximum likelihood 
estimator for \riD 2 is asymptotically unbiased, and sat- 
urates the Cramer-Rao bound on the uncertainty given 
by the Fisher information. 



A. No intrinsic dispersion 

We consider a distribution of sources with some magni- 
fication probability P{x), where x = \n\x. The apparent 
distance to the source is 
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Our job is then straightforward: we are to estimate InP/ 2 
from N independent and identically distributed values of 
In -D 2 pp j • If AT is sufficiently large (how large will be in- 
vestigated below), then we may use the Fisher informa- 
tion to determine how well we can measure In D 2 . For a 
single sample (N — 1), the Fisher information is 
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For multiple independent samples, the Fisher matrix sim- 
ply adds so that I\ n ij2 — NlP^ D2 . For large N, the un- 
certainty in In D 2 would then be 
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B. Intrinsic dispersions and measurement 
uncertainties 

We now consider the case where the error x in In D 2 __ 

a pp 

is determined not just by lensing, but also by an addi- 
tional contribution such as intrinsic dispersion (for statis- 
tical standard candles such as supernovae) or flux mea- 
surement uncertainty (which exists for any standard can- 
dle). We denote the lensing contribution by X\ and the 
additional dispersion by x-i. We assume these to be in- 
dependent with probability distributions Pi and P2, so 
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that the probability of x is given by a convolution: 



P(x) 



Pi(xi)Pi(x - xi)dxi. 



(5) 



This assumption should be true for the case where £2 
is dominated by intrinsic dispersion, since the intrinsic 
luminosity scatter of a source physically cannot depend 
on the lens alignment. It might be violated for the case 
of the measurement uncertainty since a magnified source 
will be detected at higher S/N and thus is likely to have 
a smaller fractional error on the flux; however this is 
probably only significant for the strongly lensed sources, 
which do not dominate the information integral, Eq. (J3|) . 

In this paper, the intrinsic dispersion/measurement 
uncertainty will be taken to be a lognormal distribution: 



The Fisher error, [A/T,' 1 ^] -1 / 2 , is of course always less 
than or equal to Eq. (fTU|) . We give an elementary proof 
of this inequality in Appendix [X] There we also show 
that equality holds only in the case where the magni- 
fication PDF is a T-distribution, Eq. (|A2I) . We expect 
that in practice Eq. (|4]) should be a substantial improve- 
ment over Eq. (fT0|) because the T-distribution does not 
resemble a realistic magnification PDF, since it cuts off 
exponentially at large magnifications. The T distribution 
is also far more symmetric than "real" PDFs: it always 
has a normalized skewness 
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equal to S' 3 = 2. 
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with the offset — er 2 2 /2 designed to ensure (e X2 ) = 1. 

The Fisher information for the convolved distribution, 
and for its improvement ratio, can be obtained from the 
usual formula, Eq. (|3J). 



III. NUMERICAL RESULTS 

We have now completed our survey of the theory; it 
is now time to actually evaluate the Fisher information 
for realistic magnification PDFs. We first describe the 
construction of the lensing magnification PDFs and then 
display results. Finally we simulate the effect of a finite 
number of sources on the maximum likelihood estimator. 



C. Estimators 

In the limit of large N, the maximum likelihood esti- 
mator for In D 2 achieves the Fisher information errors. 
This estimator is given by the implicit equation 
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is a weight function. 

This can be compared to the "conventional" estimator 
based on flux-averaging [8J, i.e. based on conservation of 
mean surface brightness (fi) = 1, which implies (D~ 2 p ) = 
D~ 2 . This approach gives another distance estimate, 
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where the subscript "C" is used to denote the conven- 
tional estimator. Note that Eq. © is model-independent 
in the sense that no functional form for the magnification 



PDF is assumed, 
viation of Eq. dU) 
tainty is 
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A. Lensing PDFs 

We use the stochastic universe method presented in 
Holz & Wald to calculate the lensing PDFs. This 
method calculates the full (weak and strong) lensing dis- 
tributions utilizing a Monte Carlo code: the universe in 
the vicinity of a photon path is generated randomly, and 
the lensing effects from the matter distribution are cal- 
culated analytically. We approximate the matter in the 
universe as pure dark matter smoothly distributed in 
NFW halos with the halo masses drawn from the 
Sheth-Tormen mass function [20I ] . and with cosmologi- 
cal parameters D, m = 0.28, = 0.72, h — 0.7, and 
<Tg = 0.79. The lensing distributions are relatively insen- 
sitive to both the details of the lenses and the values of 
the cosmological parameters [2l| . 

The Fisher analysis requires that the magnification 
PDF be smooth, since Monte Carlo noise results in a 
spurious, positive definite contribution to Eq. ([3]). We 
have used several versions of the smoothing procedure. 
The default procedure (used for most of the results in 
this paper unless otherwise specified) is to bin the lens- 
ing PDF in bins of width A/i = 10~ 3 . Then a triangle-hat 
smoothing kernel is used, i.e. 
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Since more smoothing is needed in the tail of the distri- 
bution than the peak, we generate two smoothed distri- 
butions Pi and P2 with different values of the smoothing 
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Magnification distributions 



Distance errors as a function of redshift and intrinsic dispersion 
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FIG. 1: The magnification distributions P(/i) at z — 0.5, 1, 
and 2. 

Si and S2. The distributions are combined according to 

P 2 + cPi 



oth 



1 + C 
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where c = e 50 (M-i-03) _ This results in an effective smooth- 
ing of Pi at large \i and P2 at small /1 with smooth in- 
terpolation. For z — 0.5 we choose (S\,S-2) = (2,15); 
for 2 = 1, (5,20); and for z = 2, (10,30). The Monte 
Carlo PDF is generated only out to /1 = 2; above this, 
we assume a P(/i) cx scaling as appropriate for large 
magnifications (near a caustic). This matters little since 
this region contributes little to Eq. ([3]). The smoothed 
distributions at z = 0.5, 1, and 2 are shown in Figure [TJ 

We have tried other methods of smoothing to ensure 
robustness. For example, we have tried re-computing the 
Fisher integral, Eq. ©, by fitting P(/x) with least-squares 
quadratic functions in intervals of width Afi = 0.02 (at 
fi > l + 0.04z) or 0.01 (at fi < l + 0.04z), and using the fit 
to analytically compute dP/dx. The integral is chopped 
off at the 0.1-percentile point of the distribution to avoid 
spurious effects (such as fits that pass through zero) since 
the quadratic polynomial is not a good approximation 
near the minimum value of fi. This procedure led to un- 
certainties that differed by at most 6% from our fiducial 
smoothing procedure. 

For completeness, we have also utilized a Savitzky- 
Golay smoothing filter, with width A/j = 0.05, and have 
found results differing by ~ 4% from our fiducial smooth- 
ing. 
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FIG. 2: The solid curves show the Fisher matrix error per 

source, [-q V^] - 1//2 j as a function of intrinsic dispersion a X2 
for source redshifts of 0.5 (bottom), 1.0, and 2.0 (top). The 
dashed curves show the error per source using the flux- 
averaging approach, Eq. {TJ. 



Eq. (JTJ). For the flux-averaging error, we have used the 
approximation er M w 0.088z [8(. The results are shown for 
3 redshifts, z — 0.5, 1.0, and 2.0, and for a range of intrin- 
sic dispersions a X2 . For large intrinsic dispersion, lensing 
adds negligible additional dispersion and the error per 
source is a X2 . For small a X2 , the lensing dispersion dom- 
inates, and we see that the Fisher matrix errors (solid 
curves) are factors of 2-3 below the flux-averaging errors 
(dashed lines). 

In the case of the flux-averaging method, the variance 
per source is as noted above given by the quadrature sum, 
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In the case of the centroiding method, no such exact re- 
lation holds. However, it turns out that the quadrature 
sum approximation 
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has an error of at most 10% over the range of redshifts 
z > 0.5 probed (the maximum error is at low redshift 
where the magnification PDF is most non-Gaussian). By 
construction, Eq. (|15p is exact when one or the other 
source of error dominates. 

For an intrinsic scatter a X2 0.15 appropriate to Type 
la supernovae, the intrinsic scatter is dominant over the 
lensing even for the flux-averaging method for z < 1.7. 
At z = 1.7 the centroiding method reduces the effective 
error per source [-fuHja] -1 ^ 2 from 0.22 to 0.16, which is 
a modest improvement but not nearly what one finds for 
low dispersion gravitational wave sources. 



Fisher results 



C. Finite sample size 



In Figure [21 we show the Fisher matrix error per 
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as well as the flux-averaging error, 



The maximum likelihood estimator is known to achieve 
the Fisher uncertainty only in the limit of large numbers 
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of sources. For a finite number of sources per redshift 
bin, Eq. (J7J may be biased and may have a larger error 
than the Fisher information would suggest. The bias in 
the estimator for In D 2 can be computed and subtracted 
using Monte Carlo simulations, assuming that P(n) is 
known. However, the uncertainty in In D 2 may be larger 
than the Fisher estimate. Here we use Monte Carlo simu- 
lations to measure the scatter in In D 2 , and show that for 
N > 4 sources the Fisher estimate is accurate to within 
< 10%. 

We have constructed our Monte Carlo simulations for 
any source redshift z and intrinsic dispersion a X2 by first 
drawing N random deviates from the P(/x) distribution. 
For each fii, we obtain an estimated luminosity distance 
D 2 = D 2 ^ 1 . Then we maximize 

N 

\nL{D) = ^P(x = \nD 2 -\nD 2 ). (16) 

i=l 

Then q = In D 2 —In D 2 has a probability distribution that 
depends only on N and P(/i) (D trivially drops out). A 
95% confidence interval on In D 2 can be obtained by tak- 
ing InD 2 and subtracting the 2.5th or 97.5th percentiles 
of the distribution of q, i.e. at 95% confidence 

InD 2 - go.975 <\nD 2 < \nt) 2 - q .Q 25 , (17) 

where q a is defined by J_° P{q)dq = a. The width of 
the confidence interval is go. 975 ~ <7o.025- 

In comparison, the usual assumption of Gaussian er- 
rors with width given by the Fisher calculation would 
suggest that the width of the 95% confidence interval 
would be 2 • 1.960[ATJ 1 ^ 3 ]- 1 /2 ) where 95% of the proba- 
bility in a unit normal distribution lies between ±1.960. 

The widths resulting from the full Monte Carlo pro- 
cedure are compared with the Fisher calculation in Fig- 
ure [3l As expected, the simulated errors are larger than 
the Fisher prediction. However the discrepancy drops 
rapidly for N > 4. This is because even 4 events are 
usually sufficient to identify and reject the strongly mag- 
nified sources. 

One would intuitively expect that the Fisher errors are 
approached more rapidly in the presence of intrinsic dis- 
persion because this results in a P(x) that is more nearly 
Gaussian. Indeed, this is what we find in simulations. 
For example, at z = 1 and N — 4, we find that with 
no intrinsic dispersion the 99% confidence region is 1.12 
times wider than the Fisher calculation suggests. This 
factor drops to 1.09 if we impose an intrinsic dispersion of 
a X2 = 0.02, to 1.07 at a X2 = 0.05, and 1.04 at a X2 = 0.1. 

IV. COSMOLOGICAL CONSTRAINTS 

We now consider the possible cosmological constraints 
from gravitational wave sources. We begin by describing 
our methodology for computing parameter constraints 
(Sec- IIV A"|) . We then turn to two specific examples: BBO 
(Sec. ITVBl and LISA fSec. lTVCl) . 
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FIG. 3: The full width of the 95% and 99% confidence 
ranges for InD 2 for source redshifts of 0.5 (bottom curve), 
1.0, and 2.0 (top curve), as a function of the number of 
sources N. Zero intrinsic dispersion is assumed. The points 
are the results from Monte Carlo simulations. The dashed 
lines are the Fisher predictions assuming Gaussian errors, i.e. 
2 ■ 1.960[NI^ D2 ]- 1/2 and 2 • 2.576[Af7 1 ( n 1) D2 ]~ 1/2 respectively. 
Note the good agreement of the Monte Carlo and Fisher re- 
sults for N > 4. This plot used 10 4 simulations. 



A. Forecasting methodology 



It is straightforward to generalize Eq. ([3]) to N sources, 
at a range of range of redshifts (z\, Z2, • • • , Zjv)j a- n d to a 
cosmological model depending on Np parameters. We 
denote the cosmological parameters by {\ a }a=i- Let 
■Capp,i represent the N measured luminosity-distances, 
and let Xi = (lnD(zj) 2 — lnD 2 pp ; ). The magnifications 
for each source should be very close to statistically in- 
dependent, since a gravitational wave detector sees the 
whole sky. (Note that this is unlike the case of a super- 
nova survey with an optical telescope, which inherently 
has a narrow field of view and hence depending on the 
survey strategy may suffer from correlated magnifications 
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Thus the ./Vp-dimensional Fisher matrix is 



I a p= I P(x 1 )P(x 2 )---P(x N ) 
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where we have used the fact that P(xi,X2,- • ■ ,x N ) = 
P(xi)P(x 2 ) ■ ■ ■ P(xn), since the N measurements are in- 
dependent. Using the fact that / P(xi)dxi = 1 and hence 
that J P(xi)[d\nP(xi)/dX a ]dxi = 0, it is easy to see in 
the double-sum over i and j, the terms with i ^ j are all 
zero. Hence our expression for I a p reduces to 
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Of course, dxi/d\ a is just d(\nD 2 )/dX a . Thus we arrive 
at our final expression for the Fisher matrix: 
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The information for a single source l[^ D 2 (zi) is simply the 
integral, Eq. ([3]), where the probability distribution for 
x contains lensing noise and (if significant) measurement 
noise as well. 

To rapidly estimate I^ D2 (z) for any z (up to z = 3), 
we (i) calculated the magnification distribution P z (fj), 
for redshifts z € {0.5, 1, 1.5, 2, 2.5, 3} using the method of 
Holz & Wald 0|, (ii) used (smoothed versions of) these 

distributions to calculate I^ D 2 (z) for these 6 redshifts, 
using Eq. ©, and then (iii) fit these results to a smooth 
curve that has the correct asymptotics (going to as 
z — >• and going to a constant as z — > oo). We find the 
following to be a suitable fit: 
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where C = 0.066, j3 = 0.25, and a = 1.8. This function 
is shown in Figure |U Note that we have not explored its 
validity at z > 3. 

At low redshifts, the peculiar velocity error dominates; 
assuming a width of 300 kms -1 (e.g. [23]) gives an ad- 
ditional contribution to a X2 of twice 300 kms -1 divided 
by the Hubble velocity cz, which is 0.002z~ 1 . We have 
approximated this error by a quadrature-sum with the 
lensing + measurement noise, which our tests of Eq. (|15[) 
suggest will not be in serious error. 

In some cases, a very large number of sources will be 
observed (possibly of order 10 5 for BBO). In such cases, 
it is appropriate to bin the sources into redshift slices 



Lensing Information versus Redshift 




FIG. 4: The uncertainty per source on the Hubble plot, 

[^l'ni 2 ( z )] _1 ^ 2 ' as a function of redshift. The points are ob- 
tained from evaluation of the information integral over P{/-i), 
while the curve is the fit of Eq. (1221) . 



as is often done for supernova forecasts [24]. Given a 
redshift distribution H(z), we bin sources into redshift 
slices of width Az = 0.1, and the number of sources in 
each bin is computed according to Ni — N SIC H(z)Az. 
We have cut off the bins at redshifts below z min , where 
z m i n is defined such that we expect 1 source at z < z mini 
i.e. J Zml ° H(z)dz — N~*. This is to prevent a mission 
that collects a small number of low-z sources from tak- 
ing advantage of a highly precise "constraint" obtained 
locally (z <C 1) from e.g. 0.01 sources with ultra-precise 
distances. 

Our parameter space {X a } includes the present-day 
densities of baryons £lhh 2 , matter fl m h 2 , and dark en- 
ergy fid c h 2 , as well as the spatial curvature . The 
equation of state (pressure: density ratio) of the dark en- 
ergy is described by the 2-parameter model 



w de (a) 



Pdc(a) 
Pde(a) 



w + w a (l - a), 



(23) 



where the parameters are (wo,w a ) and the fiducial "cos- 
mological constant" model has wq — — 1 and w a = 0. We 
also include the primordial spectrum of Gaussian adia- 
batic cosmological perturbations, assumed to be a power 
law (2 parameters: amplitude and spectral index), which 
are required when combining GW data with the CMB or 
weak lensing; and the absolute magnitude of a Type la 
supernova, required when including the supernova Hub- 
ble diagram. This gives a 9-dimensional parameter space. 
Note that the Hubble constant Hq is not an additional 
parameter since it is determined by fi m /i 2 , fl^ch 2 , and 
fl K h 2 . 

We run our forecasts both internal to the GW method, 
and in combination with other methods of probing cos- 
mology; the latter cases include the "Stage III" (i.e. near- 
term ground based) results for supernovae, weak lensing, 
and baryon-acoustic oscillations, and the Planck CMB 
constraints, as forecast by the Figure of Merit Science 
Working Group (FoMSWG) 0- 
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We compute the Figure of Merit (FoM) defined by the 
Dark Energy Task Force (DETF) [25[ for several combi- 
nations of future data sets. This FoM is defined to be 
proportional to the inverse-area of the error ellipse in the 
{w Q ,w a ) plane, i.e. 



No intrinsic dispersion 



FoMdetf 



a(w )a(w a )y / l - p 2 (w ,w a ) 



(24) 



where p(wo,w a ) is the correlation coefficient. The DETF 
Figure of Merit is not a unique (or even necessarily the 
best) way to present the performance of a dark energy 
experiment - see the discussion in the FoMSWG report 
[24| - but it is well suited to our objective here, which is 
to show that our magnification PDF centroiding method 
leads to significantly tighter dark energy constraints from 
GW sources. 



B. Example: BBO 

We consider a population of sources with the red- 
shift distribution of Ref. [5] , appropriate to neutron star- 
neutron star (NS-NS) binaries. Two cases are consid- 
ered for the error of the distance determination: an 
ideal case (IDEAL, a X2 = 0), and an error of \Az% 
(NSNS, a X2 — 0.028z), with the latter appropriate to 
BBO parameter forecasts [5|. For each of these cases, 
we consider two subcases for the distance determina- 
tion, the flux-averaging method (AVE) and the cen- 
troiding method described here (CEN). For the centroid- 
ing case, we used the method described above, while 
for the flux-averaging cases, the lnZ? 2 uncertainties per 
source are 0.088z (IDEAL. AVE) or V0.088 2 + 0.028 2 z 
(NSNS. AVE). 

In Figure El we show the DETF Figure of Merit for the 
combined constraints. We have not included any system- 
atics in the GW constraints. This is of course an opti- 
mistic case, and it is not yet clear whether the systematic 
errors can be made negligible for a BBO-class mission. 
For example, while the physics of the GW source (the 2- 
body problem in general relativity) is "clean," there are 
possible systematic errors associated with (i) the theo- 
retical predictions for the magnification PDF P(n), par- 
ticulary associated with baryonic physics on small scales; 
and (ii) the strain calibration of a BBO-type detector Q. 

We can see that an improvement of a factor of ~ 2 in 
the FoM is possible with ~ 130 well-measured sources, 
and an order of magnitude improvement (comparable to 
that promised by various Stage IV projects such as the 
more grandiose versions of JDEM |42|) is possible with 
~ 1500 sources. We also see that using the flux averaging 
rather than the centroiding results in a factor of ~ 5 
increase in the number of sources required to reach a 
given DETF FoM for the IDEAL case (and a factor of 
- 3 for the NS-NS case). 
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FIG. 5: The DETF FoM as a function of the number of grav- 
itational wave sources N STC used. We also include both the 
Planck mission and next-generation ground-based dark en- 
ergy projects (Stage III). The highest N SIC value plotted cor- 



responds to N al 
BBO. 



3 x 10 , the rough number expected for 



C. Example: LISA 

As a second application of the main ideas in this paper, 
we consider how well cosmological parameters might be 
constrained by LISA observations of coalescing massive 
black hole binaries (MBHBs). This question has been 
considered by several authors, including (ill [26h 29] . with 
the importance of weak lensing as the dominant effect in 
limiting LISA's distance-measurement accuracy first be- 
ing stressed by Hughes & Holz [30T |. The main result of 
this paper - that previous analyses considerably under- 
estimated the improvement in D^-accuracy that comes 
from combining several measurements - suggests a re- 
examination of the LISA case. 

To provide some context: LISA will be capable of 
detecting merging MBHBs with masses in the range 
- 10 3 - 1O 6 M out to z ~ 20. Estimates of MBHB 
merger rates vary by several orders of magnitude, de- 
pending mostly on the fraction of dark-matter halos that 
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have MBHs in their cores at redshifts z ~ 10 — 20. For ex- 
ample, merger-tree models due to Volonteri predict that 
LISA should detect ~ 30 MBHBs/yr, mostly at high red- 
shift (z > 4) [31(. LISA's angular resolution will typi- 
cally be a few degrees or worse, so to uniquely identify 
the host galaxy will generally require some corresponding 
electromagnetic outburst. Several possible mechanisms 
for generating outbursts have been explored, including 
(i) excitation of a shared accretion disk due to the rapid 
mass loss and/or velocity kick when the binary merges 
(a consequence of the energy and momentum carried off 
in GWs) [3214341 ] . (ii) a steep rise in the accretion rate 
starting months to years after the merger [35| , and (iii) a 
pre-merger burst due to shepherding of the disk around 
one of the progenitor black holes [36| ■ But for accurate 
knowledge of the intensity, spectrum, and time-profile 
of such electromagnetic outbursts, we may well have to 
wait until LISA flics. Lacking robust predictions regard- 
ing electromagnetic outbursts, the LISA community has 
tentatively adopted the criterion that an MBHB merger 
is promising for precise localization if the LISA lcr er- 
ror ellipse on the sky is < 10 deg 2 , which is roughly the 
field of view of the planned Large Synoptic Survey Tele- 
scope (LSST) [43[. Applying this criterion to Volonteri's 
population models, for example, one finds that ~ i of 
LISA's detected MBHBs could be positioned to within 
< 10 sq. degrees. Even if only ~ 20% of these "promis- 
ing" events could actually have their redshift determined, 
this would still lead to of order 10 points on the Dl — z 
curve where the luminosity distances follow from funda- 
mental physics (the 2-body problem in GR). Errors in Dl 
due to noise (both instrumental noise and the confusion 
background from ~ 3 x 10 8 compact galactic binaries) 
will typically be only *~ 1 — 2%, even before incorporat- 
ing the extremely accurate sky-localization provided by 
an EM counterpart. Using the precise EM localization 
will typically reduce this uncertainty by a factor ~ 2 [2^] . 
Therefore we are in a regime where WL magnifications 
strongly dominate the Dl errors. 

Given the large uncertainties, in this paper we adopt 
a very simple model for the z-distribution of localiz- 
able sources, which is as follows. We take the rate of 
MBHB mergers (per unit comoving volume, per unit 
proper time) to be some constant n, with the universe 
evolving according to a flat ACDM model, with our fidu- 
cial values (fi A = 0.744, fl m = 0.256). Then the redshift 
distribution of the of the binary sources whose GWs are 
arriving at LISA (over an observation time T Q b s ) is: 



dN r 2 {z) 
_= 47 ™T obS(i + z)H(z) , 



(25) 



where r(z) is the comoving distance to redshift z and 
H(z) is the Hubble rate. 

We restrict attention to mergers at z < z max = 3, and 
we assume that some constant fraction F of all merg- 
ers in this redshift range can be associated with an EM 
outburst that identifies the host galaxy. The rough justi- 
fication for limiting our attention to redshifts z < z max is 



that as z increases it clearly becomes harder to identify a 
counterpart, both because the GW SNR is lower (which 
increases the size of the error box) and because any EM 
outburst is fainter (and at z > 7.5 is completely ob- 
scured in the LSST bandpasses by intergalactic Lyman-a 
absorption); we crudely model this fall-off by a Heavi- 
side function. We generate parameter constraints by the 
method described in Section IIV Al 

The resulting parameter constraints are shown in Fig- 
ure [S] for both the case of 10 and 30 usable electromag- 
netic counterparts. The addition of the GW constraint 
does not significantly improve the Stage III DETF Fig- 
ure of Merit - for 30 sources it raises it from 116 to 130 
(although we note that the investigation of GW dark en- 
ergy constraints is still in its early days and further im- 
provements may be possible). However, one may assess 
the robustness of an overall dark energy program in part 
by examining how well one can do with each dark en- 
ergy technique pEj . We therefore show in Figure [6] con- 
straints for the CMB+SN+GW, CMB+WL+GW, and 
CMB+BAO+GW cases. The gravitational wave con- 
straints make large improvements when combined with 
the supernovae (they partially break the w a — degen- 
eracy by extending the Hubble diagram to higher red- 
shifts) and weak lensing. Less improvement is seen with 
the BAO model because the BAO already provide some 
distance constraints in the z > 1 range. As one can see by 
comparing the top and bottom rows of the figure, the pa- 
rameter constraint improvements are much more impres- 
sive when using the centroiding than the flux-averaging 
method. 



V. DISCUSSION 

The luminosity distance-redshift relation is one of the 
oldest and most fundamental probes of cosmology, and 
future gravitational wave detectors offer the possibility 
of measuring it accurately using binary inspirals whose 
luminosity can be calculated directly from measured pa- 
rameters and fundamental physics. These sources are 
however affected by weak gravitational lensing by inter- 
vening inhomogeneities in the cosmic mass distribution. 
This introduces changes of typically a few percent (but 
occasionally much larger) in the flux, while not signifi- 
cantly affecting the redshift, and thus provides a source 
of noise in the D(z) relation. We have shown in this pa- 
per that exploiting the full power of the likelihood func- 
tion can reduce this noise: the noise in the D{z) relation 
is not limited by the lensing dispersion divided by the 
square root of the number of sources, but rather can be 
less by a factor of 2-3 if one centroids the distribution of 
apparent distances using the known non-Gaussian form 
of the lensing magnification PDF. 

We have not discussed here the systematic errors as- 
sociated with using large numbers of gravitational wave 
sources for precision low-redshift cosmology, as suggested 
for BBO. While the signal itself is expected to be clean, 
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FIG. 6: The constraints on the (wo,w a ) model. The solid ellipses are the forecast 68% confidence contours (Ax 2 = 2.28) for 
Planck plus the indicated cosmological probe (SN, WL, or BAO) at Stage III using the FoMSWG Fisher matrices [24| . The 
dashed and dotted ellipses incorporate LISA constraints assuming either 10 usable sources at z < 3 (dashed) or 30 sources 
(dotted). The upper panels shows results using the centroiding technique (CEN), while the bottom panels use flux averaging 
(AVE). The horizontal axis is plotted as w(z = 0.5) = wo + ^w a instead of wq in order to avoid long diagonal contours; note 
that this transformation preserves the areas of the ellipses. 



there are potential sources of systematic error. Some of 
these, such as strain calibration, coherent peculiar veloc- 
ities at low redshift [37], H|| , and host redshift misidentifi- 
cation, exist irrespective of the method used to estimate 
the true D(z) from a collection of weakly lensed GW 
sources with their apparent fluxes and redshifts. How- 
ever, the issue of uncertainties in the magnification PDF 
is worth discussing here. It may seem at first glance that 
the flux-averaging technique is more robust than the cen- 
troiding technique described here, because it relies only 
on the flux conservation constraint J °° P(/i)d/i = 1. This 
may not be the case for three reasons. One is that in 
order to remove misidentified host galaxies, it is likely 
that a BBO-type mission would reject outliers from the 
D(z) relation [5]. This outlier rejection would elimi- 
nate the the tails of the magnification distribution with 
|ln/i| > (0.4 In 10) AM, where AM is the half-width of 
the cut in magnitudes. Since P(/x) is highly asymmet- 
ric, with the large-/x tail much stronger than the small- 
fi tail, it follows that outlier rejection will result in the 
conditional probability J Q P(/i|accept)<i/j < 1 and hence 
give a positive bias in D(fi). This can be corrected, but 
it requires knowledge of the contribution to J Q P(^)d^i 
from the strong-magnification tail. A second reason has 



to do with strong lensing: flux conservation implies that 
the magnification satisfying J* °° P(/i)efyi = 1 is the total 
magnification of all of the images. However, since strong- 
lens time delays are often measured in months (and even 
longer if the strongly de-magnified central image is signif- 
icant), it is likely that for many BBO sources there will be 
additional images whose time delay places them outside 
the BBO mission lifetime. Finally, in obtaining a success- 
ful host redshift (or identifying a source with the correct 
host rather than another object nearby on the sky), there 
is likely to be a bias in favor of brighter sources, which 
results in a success probability that has some dependence 
on the magnification, P(success) ~ /i* 3 . The presence of 
lensing dispersion then results in a bias in the mean flux 
of ~ /3Var/i; since Var fi is of order 10 -2 at z ~ 1, such 
biases will likely be far above the BBO statistical errors, 
and will have to be corrected using knowledge of P{p)- 
Therefore, even the flux-averaging method is sensitive to 
the particular distribution P(/i). The problem of how 
well P(/x) can be determined via theory (particularly in 
the presence of baryonic effects), or reduced to a para- 
metric model whose parameters can be simultaneously fit 
using BBO, is left to future work. 

In this paper, we have not attempted to reduce the 
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lensing dispersion by using external information, i.e. only 
the shape of P(p-) was assumed. The main external 
sources of information that are commonly considered are 
weak lensing shear maps and galaxy maps. In principle 
the way one would incorporate this information would 
be to write a conditional probability density, e.g. P(/i\g), 
where g represents the galaxy density field. Galaxy maps 
have been shown to be helpful because they contain in- 
formation about the small scales that dominate the lens- 
ing variance [H|-[l5j]; however the conditional probabil- 
ity distribution P(fj,\g) may be very hard to determine 
theoretically at the required precision. For example, de- 
spite recent advances in determining the relation between 
galaxy luminosity and host halo mass (a key quantity of 
interest if one is trying to infer the matter density field 
from galaxy observations) using clustering and lensing 
data [3 3 . |40j , a measurement of the scatter in this rela- 
tion is not yet possible, and the full distribution of this 
scatter - required if one is going to compute P(fi\g) - is 
woefully underconstrained. Nevertheless, for a mission 
such as LISA that may have only a limited number of 
usable sources and hence may be dominated by statis- 
tical errors due to weak lensing, this may be a useful 
approach. The weak lensing field is sensitive only to the 
matter distribution, and so one could imagine that it 
would be profitable to utilize the smoothed convergence 
field, K sm , and attempt to centroid the conditional den- 
sity P(n\k S m)- WL has traditionally been viewed as not 
useful for de-lensing of GW sources because most of the 
lensing variance comes from small scales where weak lens- 
ing measurements are noisy [12| . This conclusion should 
be revisited in future work using the centroiding tech- 
nique; in particular, this small-scale structure contributes 
strongly to the high- magnification tail of P(/x|re sm ), and 
it is not yet known what happens to the Fisher informa- 
tion, which depends largely on the width of the peak of 
P(fj,\k am ). 
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Appendix A: Information improvement 

This appendix is dedicated to proving the information- 
bounding inequality: 



rr(l) i-l/2 < 



(Al) 



for any distribution such that (fi) = 1. We also show 
that equality holds only for the T-distribution: 



a-l -ufl 



We begin by considering the functions 



f(x) = (e* - l)y/P@j 



and 



-2f y/P(x). 

dx 



(A2) 



(A3) 



(A4) 



We now use the Cauchy-Schwarz inequality, 
[J f(x)g(x)dx} 2 < [J f 2 (x)dx][J g 2 (x)dx}. It is readily 
apparent that a 2 ^ — J f 2 (x)dx since /j, — e x and (fi) = 1. 
We can also see that 



g 2 (x)dx 



2 Tx^ 



dx 



VP(^)— lnPO) 



dx 



L lnD 2 ' 



(A5) 



Finally, 



f(x)g(x)dx = -2 j {e x -l)y/P{x) — y/P{x)dx 



(e x - 1) — P(x)dx 
dx 



e x P(x)dx = (/i) = 1, 



(A6) 



where in the third equality we have used integration by 
parts and noted that the surface terms vanish since in 
order to be normalized the probability distribution must 
vanish faster than x _1 as x — > ±oo. The Cauchy-Schwarz 
inequality then states 



1 < crf.L 



2r(l) 



M In D 2 5 



(A7) 



thereby proving Eq. (jAll) . 

Equality holds if and only if g{x) — af(x) for some 
constant a. Examining Eqs. (|A3|) and (IA4|) shows that 
equality is thus equivalent to a first-order ordinary dif- 
ferential equation for yj P(x), 



2- V / PW=a(e I -l)Vi 5 W, (A8) 
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which has the solution 

P(x) oc exp[-a(e x - x)}. (A9) 

Using P(fi) = P(x)dx/dfi = P(x)//j,, we may write this 
as a function of fi: 

P(fi) oc n^e-^. (A10) 
This equation is easily normalized and is given by 



Eq. (IA2[) . By inspection its first moment is indeed 
(fi) = 1, and its variance is Var/x = a -1 . 

We note that for small variance (large a), the re- 
distribution (Eq. IA2[) approaches a Gaussian. This is a 
direct consequence of the central limit theorem since the 
T-distribution is simply the reduced-x 2 distribution, i.e. 
X 2 /Ndof for Adof — 2a degrees of freedom, and hence 
represents the distribution of sample averages of N^ a { 
squared unit Gaussian deviates. 
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